/*******************************************************************************

                                WEEDFO.C

  Title:        Speden weed procedure for fobs files. 
  Author:       Hanna Szoke
  Date:         12/30/2004 
  Function:     Weed out near-duplicate structure factors.
		It is assumed that the file has previously been sorted,
		at least to the extent that the h values are clumped,
		using, e.g.  'sort -n +1 -2 [in_filename] >[out_filename]',
		or even just 'sort [in_filename] >[out_filename]'.

*******************************************************************************/

#include	"util.h"
#include	"cellparams.h"
#include	"dims.h"
#include	"sp.h"

float	foround = 0.1 ;
char	weed_filename[MAXSTRING], dups_filename[MAXSTRING] ;
FILE	*fp_weed, *fp_dups ;

void
weedfo_main(int argc, char *argv[]) {

     int	HDFdatalength(char *, char *) ;
     int	readfobs(char *) ;
     int	read_HDFfobs(char *) ;
     int	weedit_fo(int) ;
     void	writeFolist(char *) ;
     void	write_HDFfobs(char *) ;
     int	checkSigmas(real *, int) ;

     char	infilename[MAXSTRING] ;
     char	outfilename[MAXSTRING] ;
     int	NoDup ;

     if (argc > optind+5) 	/* too many arguments */
	  ballpark(caller) ; 
     else if (argc > optind+3) {
          sprintf(infilename, "%s", argv[optind+2]) ;
          sprintf(outfilename, "%s", argv[optind+3]) ;
          if (argc == optind+5)
               foround  = (float) atof(argv[optind+4]) ;
     } 
     else {
	  prompt("\nPlease enter 2 filenames - input and output: ") ;
          sscanf(terminp, "%s  %s", infilename, outfilename) ;
     }

     sprintf(weed_filename, "%s.weed", infilename) ;
     sprintf(dups_filename, "%s.dups", infilename) ;

     if ((fp_weed = fopen(weed_filename,"w")) == NULL) {
          sprintf(message, "Cannot open %s", weed_filename);
          Fatal(message) ;
     }
     if ((fp_dups = fopen(dups_filename,"w")) == NULL) {
          sprintf(message, "Cannot open %s", dups_filename);
          Fatal(message) ;
     }

     hello(caller) ;

     readBasicInput() ;
     if (HDFIN)
          Nfob = HDFdatalength(infilename, "/fobs/Nob") ;
     else
          Nfob = filelength(infilename) ;

     strcpy(message, "weedfo") ;

     Folist =  (real *) e_malloc(Nfob*sizeof(real), message) ;
     Siglist = (real *) e_malloc(Nfob*sizeof(real), message) ;
     Hlist =  (float *) e_malloc(Nfob*sizeof(float), message) ;
     Klist =  (float *) e_malloc(Nfob*sizeof(float), message) ;
     Llist =  (float *) e_malloc(Nfob*sizeof(float), message) ;

     printTwice("Reading input...") ;
     if (HDFIN)
          Nfob = read_HDFfobs(infilename) ;
     else
          Nfob = readfobs(infilename) ;
     useSig = checkSigmas(Siglist, Nfob) ;

     if (!useSig)
        e_free(Nfob*sizeof(real), Siglist) ;

     printTwice("Weeding...") ;
     NoDup = weedit_fo(Nfob) ;
     Nfob = NoDup ;
     printTwice("Writing output.") ;
     if (HDFOUT)
          write_HDFfobs(outfilename) ;
     else
          writeFolist(outfilename) ;
}
int weedit_fo(N)
{

	/****************************************************************
	Check each entry against all others.  Is it a (near) duplicate 
	in (h,k,l)?  (Leaving duplicates means weighting them 
	disproportionately.)  If so, replace all such by the mean.
	Report duplicates and statistics of weeding.
	****************************************************************/

     float	roundit(float *, float) ;
     int        n, nn ;
     int	ndups, totdups = 0 ;
     float	h, k, l ;
     float	hh, kk, ll ;
     real 	this_amp, this_sig, next_amp, next_sig ;
     real	suma, sums, invss, suminvss ;

     short	*dups ;

     dups = (short *) e_malloc(N*sizeof(short), "weed_duplicates") ;

     for (n = 0; n < N; n++) 
	  *(dups+n) = 1 ;

		/* Round all indices to multiples of 'foround'	*/

     for (n = 1; n < N; n++) {
        *(Hlist+n) = roundit(Hlist+n, foround) ;
        *(Klist+n) = roundit(Klist+n, foround) ;
        *(Llist+n) = roundit(Llist+n, foround) ;
     }

	/* Identify entries that have duplicates ***/
printTwice("Ready to identify duplicates\n") ;

     for (n = 0; n < N; n++) {

        if (*(dups+n) == 1) {
           h = *(Hlist + n) ;
           k = *(Klist + n) ;
           l = *(Llist + n) ;

           for (nn = n+1; nn < N; nn++) {
   
              hh = *(Hlist + nn) ;
              kk = *(Klist + nn) ;
              ll = *(Llist + nn) ;

              if ((h == hh) && (k == kk) && (l == ll)) {
                 *(dups+n) += 1 ;
                 totdups += 1 ;
              }
              else if (hh != h) 
                 break ;
           }
        }
if ((totdups > 0) && ((n%100)*100 == n)){
     sprintf(message, "n=%d, totdups=%d\n", n, totdups) ;
     printTwice(message) ;
}
     }

	/* determine what should replace the duplicates ***/

printTwice("Ready to see how to replace near-duplicates\n") ;
     for (n = 0; n < N; n++) {

        if (*(dups+n) > 1) {
           ndups = *(dups+n) ;
           h = *(Hlist + n) ;
           k = *(Klist + n) ;
           l = *(Llist + n) ;

           this_amp = *(Folist + n) ;
           if (useSig) {
              this_sig = *(Siglist + n) ;
              sums = this_sig*this_sig ;
              invss = 1. / (this_sig*this_sig) ;
              suma = this_amp*this_amp * invss ;
              suminvss = invss ;
	      fprintf(fp_dups, 
              "INDEX   %10.2f  %10.2f  %10.2f  FOBS   %g   SIG   %g\n", 
	          h, k, l, this_amp, this_sig) ;
           }
           else {
              suma = this_amp ;
	      fprintf(fp_dups, 
              "INDEX   %10.2f  %10.2f  %10.2f  FOBS   %g\n", 
	          h, k, l, this_amp) ;
           }

sprintf(message, "Candidate %d identified \n", n) ;
printTwice(message) ;
           for (nn = n+1; nn < N; nn++) {

              hh = *(Hlist + nn) ;
              kk = *(Klist + nn) ;
              ll = *(Llist + nn) ;

              if ((h == hh) && (k == kk) && (l == ll)) {
                 next_amp = *(Folist + nn) ;
                 if (useSig) {
                    next_sig = *(Siglist + nn) ;
                    sums += next_sig*next_sig ;
                    invss = 1. / (next_sig*next_sig) ;
                    suma += next_amp*next_amp * invss ;
                    suminvss += invss ;

	            fprintf(fp_dups, 
                    "INDEX   %10.2f  %10.2f  %10.2f  FOBS   %g   SIG   %g\n", 
	  	         hh, kk, ll, next_amp, next_sig) ;
                 }
                 else {
                    suma += next_amp ;
	            fprintf(fp_dups, 
                    "INDEX   %10.2f  %10.2f  %10.2f  FOBS   %g\n", 
	                 h, k, l, next_amp) ;
                 }
              }
              else if (hh != h) 
                 break ;
           }
           if (useSig) {
              *(Folist + n) = sqrt(suma / suminvss) ;
              *(Siglist + n) = sqrt(sums / ndups) ;
              fprintf(fp_weed, 
             "%d inputs near (%g,%g,%g) replaced by amp %g, sigma %g \n", 
		ndups, h, k, l, *(Folist + n), *(Siglist + n)) ;
           }
           else {
              *(Folist + n) = suma / ndups ;
              fprintf(fp_weed, 
             "%d inputs near (%g,%g,%g) replaced by amp %g \n", 
		ndups, h, k, l, *(Folist + n)) ;
           }
         }
      }
      fclose(fp_weed) ;
      fclose(fp_dups) ;

     /*** Finally - do the weeding ***/
printTwice("Finally, weed\n") ;

     n = 0 ;
     for (nn = 0; nn < N; nn++) {
         if (*(dups+nn) == 1) {
            *(Hlist + n) = *(Hlist + nn) ;
            *(Klist + n) = *(Klist + nn) ;
            *(Llist + n) = *(Llist + nn) ;

            *(Folist + n) = *(Folist + nn) ;
            if (useSig)
               *(Siglist + n) = *(Siglist + nn) ;
            n++ ;
        }
     }

     sprintf(message, "A total of %d inputs have been removed", totdups) ;
     printTwice(message) ;
     fprintf(stdout, "For details, see .weed and .dups files\n\n") ;

     e_free(N*sizeof(short), dups) ;

     return(n) ;
}
